GPS Footprinting

2021-01-12 09:32:21

Author: Cedric Huchuan Xia (email, github)

Affiliation: Penn Lifespan Informatics and Neuroimaging Center (PennLINC)

---

This GPS Footprinting project is inspired by Finn et al. Nat Neuro (2015) and Kaufmann et al. Nat Neuro (2015). These authors studied how personal differences in functional brain connectivity can identify individuals, mature during development, alter in neuropsychiatric illness, and differ between genders. The authors referred to these individual differences as brain fingerprinting.

Here, we apply the fingerprinting technique to highly sampled GPS data in a clinical sample of youth. We are interested in how their individual mobility patterns can distinguish one another’s identity, differ between genders, and alter in psychopathological groups. In other words, can our footprint tell us apart and something about our behavior?

GPS data preprocessing was performed according to Ian Barnett’s imputation algorithm, originally published in Biostatistics (2020).

1. Setup Environment

require(ggplot2)
require(summarytools)
require(cowplot)
require(caret)
require(corrplot)
require(RColorBrewer)
require(vembedr)
require(Rmisc)
require(varian)
require(patchwork)
require(plotly)
require(Metrics)
require(dplyr)
require(ggpubr)
require(mosaic)
require(openxlsx)
require(visreg)
require(factoextra)
require(rstatix)
require(gridExtra)
require(colorspace)
require(grid)
source('~/Documents/GitHub/SmartphoneSensorPipeline/Extra/plotting_functions.R')

project_path = "~/Documents/xia_gps/"
data_path = file.path(project_path,"beiwe_output_043020")
gps_df_path = file.path(data_path,"Processed_Data/Group/feature_matrix.txt")

2. GPS Features

The current data has 41 subjects, consisting of 3317 total days, and 15 GPS features. To conserve battery life, a subject’s GPS coordinates were tracked for in a 2-min-on and 18-min-off cycle everyday using their own mobile device via the Beiwe platform. We used Ian Barnett’s algorithm to impute the missing data during the off cycles. Assuming no more data was missing due to various factors, one would generate 144 mins of GPS data per day, or 10% of total minutes in a day.

As you can see from the figure below, we can nicely reconstruct an individuals’ mobility trajectory from these data. Figure 1: A weekly view of a subject’s mobility pattern

For an even more intuitive view of the GPS data, take a look at the video here:

From these GPS traces, we extracted daily mobility features, such as max home distance, circadian routine, probability of pauses, as described in Barnett et al., Biostatistics (2020) and defined mathematically in its supplementary material.

To get a flavor of what these features look like, the first six days of GPS data for a subject are attached below.

gps_df = read.table(gps_df_path,header = T, dec = ",", )[,c(1,2,97:111)]

#Define the column types
gps_df$Date = as.Date(gps_df$Date)
gps_df[,3:dim(gps_df)[2]] = apply(gps_df[,3:dim(gps_df)[2]], 2, function(x) as.numeric(x))
head(gps_df)

NA

3. Exclude Data

The first step before further analysis is to exclude data points (days) that had excessive amount of data missing. Here is a historgram of the minutes missing for all 3317 total days.

p=minMiss_histplot(gps_df,200, "All Data")

ggplotly(p)

Figure 2: Minutes missing for all collected days. There are 3317 total days. Dashlines indicate the corresponding percentile.

3a. remove first and last days

First, we will remove the first and last days from each subject, because the application was installed during mid-day at the beginning of the study and uninstalled mid-day at the end of the study.

# loop through each subj to remove 1st and last days of gps data
gps_df_clean = data.frame() #initiate a df
for (subj in unique(gps_df$IID)){ #loop through each subj
  gps_df_subj <- subset(gps_df, IID == subj) #get gps_df per subject
  gps_df_subj <- gps_df_subj[2:(dim(gps_df_subj)[1]-1),] #remove the 1st and last days
  gps_df_clean <- rbind(gps_df_clean,gps_df_subj) #combine all subjs
}

This step removed 82 days. Now, dataset has 41 subjects, consisting of 3235 total days, and 15 GPS features.

Figure 3: Minutes missing after removing the first and last days of each subject. There are 3235 total days. Dashlines indicate the corresponding percentile.

3b. remove days at the sensitivity threshold

Next, we are removing the days with excessive data missingness. This is of course an arbitrary step. Therefore we will conduct a sensitivity analysis of the missingness threshold we set by performing the analysis across multiple different thresholds. For now, the current sensitivity cutoff is set at 1440, which amounts to only excluding those with 100% of GPS missing.

sensitivity_cutoff = 1440 # this controls the cutoff threshold

Warning messages: 1: In readChar(file, size, TRUE) : truncating string with embedded nuls 2: In readChar(file, size, TRUE) : truncating string with embedded nuls

gps_df_clean2 = subset(gps_df_clean, MinsMissing < sensitivity_cutoff)

This step removed 217 days. Now, dataset has 41 subjects, consisting of 3018 total days, and 15 GPS features.

Figure 4: Minutes missing after removing days with all data missing. There are 3018 total days. Dashlines indicate the corresponding percentile.

Figure 5: A Zoom-in plot of minutes missing distribution after removing days with all data missing. There are 3018 total days. Dashlines indicate the corresponding percentile.

subj_minmissing = gps_df_clean2 %>% group_by(IID) %>% summarise(mean = mean(MinsMissing), n = n())

subj_minmissing_plot = ggplot(subj_minmissing) + 
  geom_point(aes(x = reorder(IID, mean), y = (1440-mean)/(1440-1296))) + 
   theme_cowplot() + 
  labs(title = "Data Completeness by Subject", 
  x = "Subjects", y = "Data Completeness") + scale_y_continuous(labels = scales::percent, limits=c(0,1)) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 8), plot.title = element_text(hjust = 0.5))
ggplotly(subj_minmissing_plot)

4. Random Data Partitions

To operationalize individual footprinting prediction, we randomly split each individual’s available GPS data to two half partitions for 100 times. We used the features from the first half to determine if we can identify the same individual using their data in the second unseen half. Again, to avoid any (un)lucky random splits, we repeated data partition 100 times.

subj_days_covid = gps_df_clean2 %>% group_by(IID) %>% dplyr::tally(name = "days")$days

Error in .$dplyr::tally(name = "days") : 
  3 arguments passed to '$' which requires 2

set.seed(510)
make_subj_seq = function(gps_df_clean2,part_times) {
  subj_seq = list()
  for (subj in unique(gps_df_clean2$IID)){
    subj_data = subset(gps_df_clean2, IID==subj)
    if (dim(subj_data)[1] >5) {
    subj_seq[[subj]] <-createDataPartition(subj_data$IID,times = part_times, p =0.5)
    }
  }
  return(subj_seq)
}
part_times = 100
subj_seq = make_subj_seq(gps_df_clean2, part_times)

5. Build Correlation Matrix

Here, similar to Finn et al. Nat Neuro (2015) and Kaufmann et al. Nat Neuro (2015), we built a Pearson correlation matrix among 15 available GPS features. Below is an illustrative sample of the variables and their correlations.

# an example of subj 1, and first half
example_data = gps_df_clean2[subj_seq$`16xv6ko1`$Resample0001,3:17]
gps_cor = rquery.cormat(example_data, type = "full")

Figure 6: Correlation matrix of GPS features for a subject.

We then calculated the feature matrix for all 41 subjects in the dataset, seperately for each half data partition, and for each random split.

make_feature_matrix = function(gps_df, subj_seq, range) {
  subj_mat_1 = list()
  subj_mat_2 = list()
  for (subj in names(subj_seq)){
    print(subj)
    subj_data = subset(gps_df, IID==subj)
    subj_mat_1[[subj]] = lapply(subj_seq[[subj]], function(list) rquery.cormat(subj_data[list,range], type = "flatten", graph = F)$r)
    subj_mat_2[[subj]] = lapply(subj_seq[[subj]], function(list) rquery.cormat(subj_data[-list,range], type = "flatten", graph = F)$r)
  }
  return(list(subj_mat_1 = subj_mat_1, subj_mat_2 = subj_mat_2 ))
}

gps_clean2_feature = make_feature_matrix(gps_df_clean2, subj_seq, 3:17 )

gps_wide_matrix = data.frame()
subj_days = gps_df_clean2 %>% group_by(IID) %>% dplyr::tally(name = "Days Collected")
for (subj in arrange(subj_days,`Days Collected`)$IID){
  #if ((subj %in% subj_days$IID[which(subj_days$`Days Collected`<=10)]) == T) {
    for (part in 1:5){
      half_1 = gps_clean2_feature$subj_mat_1[[subj]][[part]]$cor
      half_2 = gps_clean2_feature$subj_mat_2[[subj]][[part]]$cor
      half_1_half_2 = c(half_1,half_2)
      gps_wide_matrix = rbind(gps_wide_matrix,half_1)
      gps_wide_matrix = rbind(gps_wide_matrix,half_2)
      }
    #}
}

gps_wide_matrix = t(gps_wide_matrix)
gps_corplot = rquery.cormat(gps_wide_matrix, type = "full",graph=FALSE)
gps_corplot$subj = arrange(subj_days,`Days Collected`)$IID

levelplot(gps_corplot$r,scales=list(draw=FALSE),col.regions = rev(rainbow(1000))[-c(1:20)], region =T, ylab.right = "Pearson correlation", main=list(label='GPS Feature Similarity'),xlab="",ylab="")

6. Match Target

Next we tried to match each target, defined by one subject’s feature matrix in the first half, to a feature matrix in the second half in the same random split. If the same individual’s feature matrix in the second half had the highest correlation among all subjects, then we assigned the match result 1, indicating a succussful match, otherwise 0, indicating an unsccussful match.

# creating a "database" against which target is to be matched with
subj_mat_1 = gps_clean2_feature$subj_mat_1
subj_mat_2 = gps_clean2_feature$subj_mat_2

calc_match_cor = function(subj_mat_1,subj_mat_2) {
  part_times = length(subj_mat_1[[1]])
  database = list() 
  for (time in 1:part_times) {
      database[[time]] = lapply(subj_mat_2, function(subjmat) subjmat[[time]]$cor)
  }
  
    # match target to database
    match_cor = list()
    for (subj1 in names(subj_mat_1)){ #loop through each subj
      # create a list of target across partitions
      target_list = lapply(subj_mat_1[[subj1]], function(part) part$cor)
      # create a match list
      for (time in 1:part_times){
        target_subj_time = target_list[[time]] #loop through each partition
        for (subj2 in names(subj_mat_2)){ #loop everyone in 2nd half
          data_subj_time = subj_mat_2[[subj2]][[time]]$cor
          match_cor[[subj1]][[as.character(time)]][[subj2]] = cor(target_subj_time,data_subj_time,use = "na.or.complete")
        }
      }
    }
  return(match_cor)
}

match_cor = calc_match_cor(subj_mat_1,subj_mat_2)

By tallying up all the successful and unsuccessful matchs in each of the 100 random splits, we calculated a distribution of match accuracy.

calc_acc_time=function(match_cor,method = "max") {
  acc_time = array()
  part_times = length(match_cor[[1]])
  for (time in 1:part_times){
    acc_time[time] = 0
    for (subj in names(subj_mat_1)){
      if (method == "max"){
        position = which.max(unlist(match_cor[[subj]][[as.character(time)]]))
      } 
      else if (method == "min"){
        position = which.min(unlist(match_cor[[subj]][[as.character(time)]]))
      }
      
      predicted_subj = names(subj_mat_1)[position]
      if (predicted_subj == subj) {
        acc_time[time] = acc_time[time] + 1
      }
    }
  }
  acc_time = acc_time/length(names(subj_mat_1))
  return(acc_time)
}
acc_time = calc_acc_time(match_cor, "max")

Over the 100 random splits, the mean match accuracy was 52.93%, with a standard deviation of 5.93%, high of 68.29%, and low of 36.59%, (95%CI: 0.52-0.54).

p_time_cor = hist_chx(acc_time, bins = 16, title = paste("Covariance Features: \n Average Accuracy Across",part_times,"Data Partitions"), xaxis = "Prediction Accuracy", yaxis = "Count")
ggplotly(p_time_cor)

Figure 7: Match accuracy distributon across 100 data splits. Mean match accuracy was 0.53 (95%CI: 0.52-0.54)

calc_acc_subj = function(match_cor, method = "max"){
  acc_subj = array()
  part_times = length(match_cor[[1]])
  for (subj in names(subj_mat_1)){
    acc_subj[subj] = 0
    for (time in 1:part_times){
      if (method == "max"){
        position = which.max(unlist(match_cor[[subj]][[as.character(time)]]))
      } 
      else if (method == "min"){
        position = which.min(unlist(match_cor[[subj]][[as.character(time)]]))
      }
      predicted_subj = names(subj_mat_1)[position]
      if (predicted_subj == subj) {
        acc_subj[subj] = acc_subj[subj] + 1
      }
    }
  }
  acc_subj = acc_subj/part_times
  acc_subj = acc_subj[-1]
  return(acc_subj)
}
acc_subj = calc_acc_subj(match_cor)

In addition to the match accuracy across subjects, we also calculated the match accuracy for each individual. Across the 41 subjects, the mean match accuracy was 37.42%, with a standard deviation of 26.71%, high of 95.8%, and low of 0%.

subj_scatter  = function(gps_df, acc_subj_vector, method){
  subj_df <- data.frame(x=names(subj_seq))
  subj_df$y = acc_subj_vector
  subj_df = subj_df[order(subj_df$y),]
  subj_acc_plot = ggplot(subj_df) + 
  geom_point(aes(x = reorder(x, y), y = y)) + 
   theme_cowplot() + 
    labs(title = paste(method,"Features \n Subject Level Accuracy Across",part_times,"Data Partitions"), 
       x = "Subjects", y = "Prediction Accuracy") +
   theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 8), plot.title = element_text(hjust = 0.5))
  return(subj_acc_plot)
}

subj_acc_plot = subj_scatter(gps_df_clean2,acc_subj, "Correlated")

ggplotly(subj_scatter(gps_df_clean2,acc_subj, "Correlated"))

Figure 8: Match accuracy distributon across 41 subjects. Mean match accuracy was 37.42%, with a standard deviation of 26.71%, high of 95.8%, and low of 0%. This shows dramatic differences in individual footprinting distinctiveness.

7. Permutation Test

To assess the statistical significance of the results above, we conducted a non-parametric permutation test. To do this, we randomly scrambled pair-wise subject-to-day linkage. We repeated this process 1000 times, and followed the same procedure as above to obtain a null distribution of average match accuracy across sample, and for each individual.

gps_df_perm = gps_df_clean2
perm_time = 1000
perm_acc_time = list()
perm_acc_subj = list()
for (i in 1:perm_time) {
  perm_part_times = 1
  print(paste("processing ...", i,"..."))
  gps_df_perm$IID = sample(gps_df_perm$IID)
  perm_subj_seq = make_subj_seq(gps_df_perm, part_times = perm_part_times)
  perm_gps = make_feature_matrix(gps_df_perm,perm_subj_seq,3:17)
  perm_mat_1 = perm_gps$subj_mat_1
  perm_mat_2 = perm_gps$subj_mat_2
  perm_match_cor = calc_match_cor(perm_mat_1,perm_mat_2)
  perm_acc_time[[i]] = calc_acc_time(perm_match_cor)
  perm_acc_subj[[i]] = calc_acc_subj(perm_match_cor)
}

Over the 1000 permutations, the mean match accuracy was 2.53%, with a standard deviation of 2.51%, high of 400%, and low of 2.5341463%.(95%CI: 0.02-0.03)

perm_acc_time_all = unlist(perm_acc_time)
q_time_cor = hist_chx(perm_acc_time_all, bins = 8, title = paste("Average Accuracy Across",length(perm_acc_time_all),"Permutations"), xaxis = "Prediction Accuracy", yaxis = "Count")
ggplotly(q_time_cor)

Figure 9: Match accuracy distributon across 1000 permutations. Mean match accuracy in permutation was 0.03 (95%CI: 0.02-0.03). This null distribution does not overlapp with the distributon using real data (Fig.7) at all.

We also calculated null distribution of match accuracy for each subject. Over the 1000 permutations, the mean match accuracy was 2.7%, with a standard deviation of 0.81%, high of 4.9%, and low of 1.1%.

perm_subj = list()
for (subj in subj_acc_plot$data$x) {
  perm_subj$val[[subj]] = sapply(perm_acc_subj, function(perm) perm[which(names(perm) == subj)])
  perm_subj$hist[[subj]] = hist_chx(perm_subj[[subj]]$val, bins = 8, title = paste(subj,": Accuracy Across \n",perm_time,"Permutations"), xaxis = "Prediction Accuracy", yaxis = "Count")
  perm_subj$acc[[subj]] = sum(perm_subj$val[[subj]])/perm_time
}

subj_df <- data.frame(x=names(subj_seq))
subj_df$y = acc_subj
subj_df$y_perm = unlist(perm_subj$acc)
subj_df = subj_df[order(subj_df$y),]
p_subj_cor_perm = ggplot(subj_df, aes(x = reorder(x, y), y = value)) + 
  geom_point(aes(y = y, col = "subject data")) + 
  geom_point(aes(y = y_perm, col = "permutation")) +
  theme_cowplot() + 
  labs(title = paste("Cor Features \n Subject Level Accuracy Across",part_times,"Data Partitions"), 
       x = "Subjects", y = "Prediction Accuracy") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 8), plot.title = element_text(hjust = 0.5))
ggplotly(p_subj_cor_perm)

Figure 8: Match accuracy null distributon across 41 subjects. Mean match accuracy in permutation was 2.7%, with a standard deviation of 0.81%, high of 4.9%, and low of 1.1%. While individuals exhibited marked differences in footprinting distinctivenss, the subject with the lowest prediction accuracy was still statistically significant against the permutation test (i.e. for subject 3e1xnphl: data: 0 % vs. permutation: 1.1%).

8. Mean Features

We calculated the mean of each GPS feature as a metric of feature variability. Similar to above, we calculated mean separately for each data partition and for each individual.

mean_gps = function(gps_df){
  mean_list = sapply(3:17,function(i) mean(gps_df[,i],na.rm=T))
  names(mean_list) = colnames(gps_df)[3:17]
  return(mean_list)
}

make_feature_mean = function(gps_df, subj_seq) {
  subj_mat_1 = list()
  subj_mat_2 = list()
  for (subj in unique(gps_df$IID)){
    subj_data = subset(gps_df_clean2, IID==subj)
    subj_mat_1[[subj]] = lapply(subj_seq[[subj]], function(list) mean_gps(subj_data[list,]))
    subj_mat_2[[subj]] = lapply(subj_seq[[subj]], function(list) mean_gps(subj_data[-list,]))
  }
  return(list(subj_mat_1 = subj_mat_1, subj_mat_2 = subj_mat_2 ))
}

gps_clean2_mean_feat = make_feature_mean(gps_df_clean2, subj_seq)

Figure 9: Feature mean across all subjects. The histogram are for data in the first random half split.

9. Variability Features

We calculated the root mean square of successive differences or RMSSD of each GPS feature as a metric of feature variability. Similar to above, we calculated RMSSD separately for each data partition and for each individual.

rmssd_gps = function(gps_df){
  rmssd_list = sapply(3:17,function(i) rmssd_id(gps_df[,i], gps_df$IID,long=F))
  names(rmssd_list) = colnames(gps_df)[3:17]
  return(rmssd_list)
}

make_feature_rmssd = function(gps_df, subj_seq) {
  subj_mat_1 = list()
  subj_mat_2 = list()
  for (subj in unique(gps_df$IID)){
    subj_data = subset(gps_df_clean2, IID==subj)
    subj_mat_1[[subj]] = lapply(subj_seq[[subj]], function(list) rmssd_gps(subj_data[list,]))
    subj_mat_2[[subj]] = lapply(subj_seq[[subj]], function(list) rmssd_gps(subj_data[-list,]))
  }
  return(list(subj_mat_1 = subj_mat_1, subj_mat_2 = subj_mat_2 ))
}

gps_clean2_rmssd_feat = make_feature_rmssd(gps_df_clean2, subj_seq)

Figure 10: Feature variability for all subjects. Variability is measured by root mean square of successive differences (RMSSD). The histograms are for data in the first random half split.

10. Predict with New Features

calc_match_vector = function(subj_mat_1,subj_mat_2,method) {
  part_times = length(subj_mat_1[[1]])
  database = list() 
  for (time in 1:part_times) {
      database[[time]] = lapply(subj_mat_2, function(subjmat) subjmat[[time]])
  }
  
    # match target to database
    match_cor = list()
    for (subj1 in names(subj_mat_1)){ #loop through each subj
      # create a list of target across partitions
      target_list = lapply(subj_mat_1[[subj1]], function(part) part)
      # create a match list
      for (time in 1:part_times){
        target_subj_time = target_list[[time]] #loop through each partition
        for (subj2 in names(subj_mat_2)){ #loop everyone in 2nd half
          data_subj_time = subj_mat_2[[subj2]][[time]]
          if (method == "cor") {
          match_cor[[subj1]][[as.character(time)]][[subj2]] = cor(target_subj_time,data_subj_time,use = "na.or.complete")
          } 
          else if (method == "rmse") {
            match_cor[[subj1]][[as.character(time)]][[subj2]] = rmse(target_subj_time,data_subj_time)
          }
        }
      }
    }
  return(match_cor)
}

match_mean_rmse = calc_match_vector(gps_clean2_mean_feat$subj_mat_1,gps_clean2_mean_feat$subj_mat_2,"rmse")
acc_time_mean_rmse = calc_acc_time(match_mean_rmse,"min")
acc_subj_mean_rmse = calc_acc_subj(match_mean_rmse,"min")

match_mean_cor = calc_match_vector(gps_clean2_mean_feat$subj_mat_1,gps_clean2_mean_feat$subj_mat_2,"cor")
acc_time_mean_cor = calc_acc_time(match_mean_cor,"max")
acc_subj_mean_cor = calc_acc_subj(match_mean_cor,"max")

match_rmssd_rmse = calc_match_vector(gps_clean2_rmssd_feat$subj_mat_1,gps_clean2_rmssd_feat$subj_mat_2, "rmse")
acc_time_rmssd_rmse = calc_acc_time(match_rmssd_rmse, "min")
acc_subj_rmssd_rmse = calc_acc_subj(match_rmssd_rmse, "min")


match_rmssd_cor = calc_match_vector(gps_clean2_rmssd_feat$subj_mat_1,gps_clean2_rmssd_feat$subj_mat_2, "cor")
acc_time_rmssd_cor = calc_acc_time(match_rmssd_cor, "max")
acc_subj_rmssd_cor = calc_acc_subj(match_rmssd_cor, "max")

11. Permutation Tests with New Features

11a. Average Accuracy

perm_vector = function(gps_df, perm_time, method_feature, method_match){

Warning message:
In readChar(file, size, TRUE) : truncating string with embedded nuls

  gps_df_perm = gps_df
  perm_acc_time = list()
  perm_acc_subj = list()
  for (i in 1:perm_time) {
    perm_part_times = 1
    print(paste("processing ...", i,"..."))
    gps_df_perm$IID = sample(gps_df_perm$IID)
    perm_subj_seq = make_subj_seq(gps_df_perm,part_times = perm_part_times)
    if (method_feature == "rmssd") {
      perm_gps = make_feature_rmssd(gps_df_perm,perm_subj_seq)
    } else if (method_feature == "mean"){
      perm_gps = make_feature_mean(gps_df_perm,perm_subj_seq)
    } else if (method_feature == "combined"){
      perm_gps_rmssd = make_feature_rmssd(gps_df_perm,perm_subj_seq)
      perm_gps_mean = make_feature_mean(gps_df_perm,perm_subj_seq)
      perm_gps_cor = make_feature_matrix(gps_df_perm,perm_subj_seq,3:17)
      subj_mat_1 = list()
      subj_mat_2 = list()
      for (subj in names(perm_gps_cor$subj_mat_1)){
          #print(subj)
          subj_mat_1[[subj]]$Resample1 = as.numeric( c(value(perm_gps_rmssd$subj_mat_1[[subj]][[1]]),value(perm_gps_mean$subj_mat_1[[subj]][[1]]),perm_gps_cor$subj_mat_1[[subj]][[1]]$cor))
          subj_mat_2[[subj]]$Resample1 = as.numeric( c(value(perm_gps_rmssd$subj_mat_2[[subj]][[1]]),value(perm_gps_mean$subj_mat_2[[subj]][[1]]),perm_gps_cor$subj_mat_2[[subj]][[1]]$cor))
      }
      perm_gps = list(subj_mat_1 = subj_mat_1, subj_mat_2 = subj_mat_2)
    }
    perm_mat_1 = perm_gps$subj_mat_1
    perm_mat_2 = perm_gps$subj_mat_2
    if (method_match == "cor") {
      perm_match = calc_match_vector(perm_mat_1,perm_mat_2, "cor")
      perm_acc_time[[i]] = calc_acc_time(perm_match, "max")
      perm_acc_subj[[i]] = calc_acc_subj(perm_match, "max")
    }
    else if (method_match == "rmse"){
      perm_match = calc_match_vector(perm_mat_1,perm_mat_2, "rmse")
      perm_acc_time[[i]] = calc_acc_time(perm_match,"min")
      perm_acc_subj[[i]] = calc_acc_subj(perm_match,"min")
    }
  }
  return(list(perm_acc_time= perm_acc_time,perm_acc_subj = perm_acc_subj))
}

p_time_mean_cor + q_time_mean_cor

Figure 11: Prediction accuracy using mean features and matched by max pearson correlation.

p_time_mean_rmse + q_time_mean_rmse

Figure 11: Prediction accuracy using mean features and matched by min root mean squared.

p_time_rmssd_cor + q_time_rmssd_cor

Figure 12: Prediction accuracy using RMSSD features and matched by max pearson correlation.

p_time_rmssd_rmse + q_time_rmssd_rmse

Figure 13: Prediction accuracy using RMSSD features and matched by min root mean squared.

combine_perm_subj = function(acc_subj_plot, perm_acc_subj) {
  perm_subj = list()
  for (subj in acc_subj_plot$data$x) {
    perm_subj$val[[subj]] = sapply(perm_acc_subj, function(perm) perm[which(names(perm) == subj)])
    perm_subj$hist[[subj]] = hist_chx(perm_subj[[subj]]$val, bins = 8, title = paste(subj,": Accuracy Across \n",perm_time,"Permutations"), xaxis = "Prediction Accuracy", yaxis = "Count")
    perm_subj$acc[[subj]] = sum(perm_subj$val[[subj]])/perm_time
  }
  return(perm_subj)
}


subj_scatter_perm = function(gps_df,acc_subj,perm_subj, method){
    subj_df <- data.frame(x=unique(gps_df$IID))
  subj_df$y = acc_subj
  subj_df$y_perm = unlist(perm_subj$acc)
  subj_df = subj_df[order(subj_df$y),]
  p = ggplot(subj_df, aes(x = reorder(x, y), y = value)) + 
    geom_point(aes(y = y, col = "subject data")) + 
    geom_point(aes(y = y_perm, col = "permutation")) +
    theme_cowplot() + 
    labs(title = paste(method,"Features \n Subject Level Accuracy Across",part_times,"Data Partitions"), 
         x = "Subjects", y = "Prediction Accuracy") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 8), plot.title = element_text(hjust = 0.5))
  return(p)
}

p_subj_mean = subj_scatter(gps_df_clean2, acc_subj_mean_cor, "Mean")
perm_subj_mean = combine_perm_subj(p_subj_mean,perm_mean_cor$perm_acc_subj)
p_subj_mean_perm = subj_scatter_perm(gps_df_clean2, acc_subj_mean_cor, perm_subj_mean, "Mean")
ggplotly(p_subj_mean_perm)

Figure 14: Subject level prediction accuracy using mean features and matched by max pearson correlation.

p_subj_rmssd = subj_scatter(gps_df_clean2, acc_subj_rmssd_cor, "RMSSD")
perm_subj_rmssd = combine_perm_subj(p_subj_rmssd,perm_rmssd_cor$perm_acc_subj)
p_subj_rmssd_perm = subj_scatter_perm(gps_df_clean2, acc_subj_rmssd_cor, perm_subj_rmssd, "RMSSD")
ggplotly(p_subj_rmssd_perm)

Figure 15: Subject level prediction accuracy using RMSSD features and matched by max pearson correlation.

12. Compare Features

p_time_cor + p_time_mean + p_time_rmssd
q_time_cor + q_time_mean + q_time_rmssd

Figure 16: Prediction accuracy across three feature sets.

p_subj_cor_perm / p_subj_mean_perm / p_subj_rmssd_perm

Figure 17: Subject level prediction accuracy across three feature sets.

13. Confounding variables

conf_scatter_plot = function(conf_var, acc_subj, col){
  if (!identical(conf_var$IID,names(acc_subj))) {
    stop("subj names do not match")
  }
  conf_df = data.frame(conf = unlist(conf_var[,col]), acc_subj = acc_subj)
  conf_scatter = ggscatter(conf_df, y = "acc_subj", x = "conf",
   add = "reg.line",  # Add regressin line
   add.params = list(color = "black", fill = "lightgray"), # Customize reg. line
   conf.int = TRUE # Add confidence interval
   ) + stat_cor(method = "pearson") +
    theme_cowplot() + theme(plot.title = element_text(hjust = 0.5)) +
    labs(title = paste(col), y = "GPS Footprint \n Prediction Accuracy", x = col)
  return(conf_scatter)
  }

subj_minmissing = gps_df_clean2 %>% group_by(IID) %>% summarise_if(is.numeric, mean, na.rm = TRUE)
sp_plots = list()
for (gps_ft in colnames(subj_minmissing)[-1]){
  sp_plots[[gps_ft]] = conf_scatter_plot(subj_minmissing,acc_subj,gps_ft)
}

sp_plots_pvals = sapply(sp_plots,function(plot) cor.test(plot$data[,1],sp_plots$Hometime$data[,2])$p.value)
sp_plots_pvals_fdr = p.adjust(sp_plots_pvals,method = "fdr")
Reduce(`+`, sp_plots)

Figure 18: Relationships between subject level prediction accuracy and individual features. Notably, there was no relationship between individual data missingness and subject level prediction accuracy.

  acc_subj_perm = acc_subj_perm[sample(length(acc_subj))]

Error: object 'acc_subj_perm' not found

Figure 19: Relationships between subject level prediction accuracy and data quantity.

14. Combine Features

gps_clean2_combined_feat = list()
for (subj in names(subj_mat_1)){
  for (time in 1:part_times){
    cor_feat1 = gps_clean2_feature$subj_mat_1[[subj]][[time]]$cor
    mean_feat1 = gps_clean2_mean_feat$subj_mat_1[[subj]][[time]]
    rmssd_feat1 = gps_clean2_rmssd_feat$subj_mat_1[[subj]][[time]]
    
    cor_feat2 = gps_clean2_feature$subj_mat_2[[subj]][[time]]$cor
    mean_feat2 = gps_clean2_mean_feat$subj_mat_2[[subj]][[time]]
    rmssd_feat2 = gps_clean2_rmssd_feat$subj_mat_2[[subj]][[time]]
    
    gps_clean2_combined_feat$subj_mat_1[[subj]][[time]] = c(cor_feat1, mean_feat1, rmssd_feat1)
    gps_clean2_combined_feat$subj_mat_2[[subj]][[time]] = c(cor_feat2, mean_feat2, rmssd_feat2)
  }
}

match_combined_feat = calc_match_vector(gps_clean2_combined_feat$subj_mat_1, gps_clean2_combined_feat$subj_mat_2, "cor")
acc_time_cb  = calc_acc_time(match_combined_feat, "max")
acc_subj_cb  = calc_acc_subj(match_combined_feat, "max")
p = hist_chx(acc_time_cb, bins = 8, title = paste("Combined Features \n by data partition"), xaxis = "Prediction Accuracy", yaxis = "Count")
q = hist_chx(acc_subj_cb, bins = 8, title = paste("Combined Features \n by subject"), xaxis = "Prediction Accuracy", yaxis = "Count")
p

Figure 20: Prediction accuracy when all three sets of features were combined.

15. Associations with Psychopathology

ids = read.xlsx(file.path(project_path,"data/clinical_data/subjecttracker_4.xlsx"))[1:length(unique(gps_df_clean2$IID)),]
#psych_score = read.csv(file.path(project_path,"data/self_report_scored_20200128.csv"))
# psych_pro = psych %>% filter(ari_proband_complete == 2)
psych_item = read.csv(file.path(project_path,"data/clinical_data/self_report_itemwise.csv"))

psych_beiwe = inner_join(ids,psych_item, by = c("BBLID" = "bblid"))

acc_subj_df = data.frame(beiweID = names(acc_subj), acc = acc_subj)
psych_beiwe_acc_subj = inner_join(acc_subj_df,psych_beiwe, by = "beiweID")


psych_sum = psych_beiwe_acc_subj %>%
            mutate(sum_als = rowSums(.[263:280], na.rm = T)) %>%
            mutate(sum_ari = rowSums(.[grep("^ari_[0-9]$",colnames(psych_beiwe_acc_subj))], na.rm = T))
psych_item

psych_sum = inner_join(psych_sum, subj_days, by = c("beiweID" = "IID"))
psych_sum$days = psych_sum$`Days Collected`
psych_sum = inner_join(subj_minmissing, psych_sum, by = c( "IID" = "beiweID" ))
psych_als = data.frame(als = psych_sum[,263:280], BBLID = psych_sum$BBLID)
psych_ari = data.frame(ari = psych_sum[,23:29], BBLID = psych_sum$BBLID)
  
 # show results

als.pca <- prcomp(na.omit(psych_sum[,263:280]), scale = TRUE)

Error in prcomp.default(na.omit(psych_sum[, 263:280]), scale = TRUE) : 
  cannot rescale a constant/zero column to unit variance

als_ari_sum = lm(acc ~ sum_als + sum_ari + days , data = psych_sum)

#als_fit = lm( acc ~ scale(pca_als) + scale(days) + scale(admin_age) + scale(admin_sex) , data = psych_sum)
#summary(lm( acc ~ scale(als_AD) + scale(days) + scale(admin_age) + scale(admin_sex) , data = psych_sum))$coefficients[2,4]

#ari_fit = lm( acc ~ scale(pca_ari) + scale(days) + scale(admin_age) + scale(admin_sex) , data = psych_sum)
#ari_als_fit = lm( acc ~ scale(pca_als) + scale(pca_ari) + scale(days) + scale(admin_age) + scale(admin_sex) , data = psych_sum)

days_plot = visreg(als_ari_sum, "days", gg = T, line=list(col="black")) + 
            ylab("GPS Footprint \n Prediction Accuracy") + xlab("Days Collected") +
            theme_cowplot() + stat_cor(method = "pearson")
sum_als_plot = visreg(als_ari_sum, "sum_als", gg = T,  line=list(col="black")) + 
            ylab("GPS Footprint \n Prediction Accuracy") + xlab("Affective Lability Score (ALS)") +
            theme_cowplot() + stat_cor(method = "pearson")
sum_ari_plot = visreg(als_ari_sum, "sum_ari", gg = T, line=list(col="black")) + 
            ylab("GPS Footprint \n Prediction Accuracy") + xlab("Affective Reactivity Index (ARI)") +
            theme_cowplot() + stat_cor(method = "pearson")
days_plot / sum_als_plot / sum_ari_plot

  psych_sum_perm$als_perm = psych_sum$sum_als[sample(length(psych_sum$sum_als))]

Error in psych_sum_perm$als_perm = psych_sum$sum_als[sample(length(psych_sum$sum_als))] : 
  object 'psych_sum_perm' not found

mod <- lm(acc~sum_als*sum_ari + days , data=psych_sum)  # just raw scores, no residuals
als_int_plot <- visreg(mod,xvar="sum_als",by="sum_ari",overlay=TRUE,strip.names=TRUE)
visreg(mod,xvar="sum_ari",by="sum_als",overlay=TRUE)

Figure 21: Associations between irritability and prediction accuracy.

16. Associations with Cognition

cnb_score = read.csv(file.path(project_path,"data/grmpy_bifactor_cnb_cleaned.csv"))
cor_traits = read.csv(file.path(project_path,"data/grmpy_corrtraits_cnb_cleaned.csv"))
cnb_beiwe = inner_join(ids,cnb_score, by = c("BBLID" = "bblid"))

17.Similarity Matrix Correlation

gps_wide_matrix_all = data.frame()

subj_seq_mat = arrange(subj_days,`Days Collected`)$IID

for (subj in subj_seq_mat){
  #if ((subj %in% subj_days$IID[which(subj_days$`Days Collected`<=10)]) == T) {
    for (part in 1:1){
      half_1 = gps_clean2_feature$subj_mat_1[[subj]][[part]]$cor
      half_2 = gps_clean2_feature$subj_mat_2[[subj]][[part]]$cor
      half_1_half_2 = c(half_1,half_2)
      gps_wide_matrix_all = rbind(gps_wide_matrix_all,half_1_half_2)
      }
    #}
}

gps_wide_matrix_all = t(gps_wide_matrix_all)
gps_corplot_all = rquery.cormat(gps_wide_matrix_all, type = "full",graph=FALSE)
gps_corplot_all$subj = arrange(subj_days,`Days Collected`)$IID

levelplot(gps_corplot_all$r,scales=list(draw=FALSE),col.regions = rev(rainbow(1000))[-c(1:20)], region =T, ylab.right = "Pearson correlation", main=list(label='GPS Feature Similarity'),xlab="Subjects",ylab="Subjects")

n_subj = length(subj_seq_mat)
gps_days_sim_matrix = matrix(NA, n_subj,n_subj)
for (i in 1:n_subj){
  for (j in 1:n_subj){
    subj_i_days = subj_days$`Days Collected`[which(subj_days$IID == subj_seq_mat[i])]
    subj_j_days = subj_days$`Days Collected`[which(subj_days$IID == subj_seq_mat[j])]
    gps_days_sim_matrix[i,j] = abs(subj_i_days-subj_j_days)/max(subj_days$`Days Collected`)
  }
}
levelplot(gps_days_sim_matrix,scales=list(draw=FALSE),col.regions = rev(rainbow(1000))[-c(1:20)], region =T, ylab.right = "Differences in Days Collected", main=list(label='Days Collected Similarity'),xlab="Subjects",ylab="Subjects")

n_subj = length(subj_seq_mat)
gps_mins_sim_matrix = matrix(NA, n_subj,n_subj)
for (i in 1:n_subj){
  for (j in 1:n_subj){
    subj_i_mins = subj_minmissing$MinsMissing[which(subj_minmissing$IID == subj_seq_mat[i])]
    subj_j_mins = subj_minmissing$MinsMissing[which(subj_minmissing$IID == subj_seq_mat[j])]
    gps_mins_sim_matrix[i,j] = abs(subj_i_mins-subj_j_mins)/max(subj_minmissing$MinsMissing)
  }
}
levelplot(gps_mins_sim_matrix,scales=list(draw=FALSE),col.regions = rev(rainbow(1000))[-c(1:20)], region =T, ylab.right = "Differences in Days Collected", main=list(label='Mins Missing Similarity'),xlab="Subjects",ylab="Subjects")

gps_sim_df = data.frame(gps = gps_corplot_all$r[lower.tri(gps_corplot_all$r)], days = gps_days_sim_matrix[lower.tri(gps_days_sim_matrix)], mins = gps_mins_sim_matrix[lower.tri(gps_mins_sim_matrix)])

days_sim_plot = ggscatter(gps_sim_df, y = "gps", x = "days",
   add = "reg.line",  # Add regressin line
   add.params = list(color = "black", fill = "lightgray"), # Customize reg. line
   conf.int = TRUE # Add confidence interval
   ) + stat_cor(method = "pearson") +
    theme_cowplot() + theme(plot.title = element_text(hjust = 0.5)) +
    labs(y = "GPS Similarity", x = "Days Collected Similarity", title = "Similarity Correlation")

mins_sim_plot = ggscatter(gps_sim_df, y = "gps", x = "mins",
   add = "reg.line",  # Add regressin line
   add.params = list(color = "black", fill = "lightgray"), # Customize reg. line
   conf.int = TRUE # Add confidence interval
   ) + stat_cor(method = "pearson") +
    theme_cowplot() + theme(plot.title = element_text(hjust = 0.5)) +
    labs(y = "GPS Similarity", x = "Mins Missing Similarity", title = "Similarity Correlation")

days_sim_plot / mins_sim_plot

18.GPS Score

gps_mean = gps_df_clean2 %>% group_by(IID) %>% summarise_if(is.numeric, mean, na.rm = TRUE)
res.pca <- prcomp(gps_mean[,2:16], scale = TRUE)
fviz_eig(res.pca)
gps_score = res.pca$x[,'PC1']
gps_score_df = data.frame(beiweID = gps_mean$IID, gps_score = gps_score)
fc_com_mean_1 = as.data.frame(sapply(fc_com_1, function(com) sapply(com, function(subj) mean(subj$V1))))
fc_com_mean_1$BBLID = as.numeric(rownames(fc_com_mean_1))
fc_com_mean_1 = inner_join(ids[,c(1,3)],fc_com_mean_1, by = "BBLID" )
fc_com_mean_1_gps = inner_join(fc_com_mean_1,gps_score_df, by = "beiweID")


fviz_pca_ind(res.pca,
             col.ind = "cos2", # Color by the quality of representation
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE     # Avoid text overlapping
             )

fviz_pca_var(res.pca,
             col.var = "contrib", # Color by contributions to the PC
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE     # Avoid text overlapping
             )

fviz_pca_biplot(res.pca, repel = TRUE,
                col.var = "#2E9FDF", # Variables color
                col.ind = "#696969"  # Individuals color
                )

19.Feature Leison

feat_names = names(gps_df_clean2[,3:17])
subj_seq_100 = make_subj_seq(gps_df_clean2, 100)
gps_ft_leison = list()
gps_df_perm = gps_df_clean2
perm_time_cut = 100
perm_acc_time_cut = list()
perm_acc_subj_cut = list()
for (i in 1:15){
  ft_rmv = feat_names[i]
  print(ft_rmv)
  range = c(3:17)[-i]
  print(range)
  # original data
  feature_nd = make_feature_matrix(gps_df_clean2,subj_seq_100,range)
  match_cor_nd = calc_match_cor(feature_nd$subj_mat_1,feature_nd$subj_mat_2)
  gps_ft_leison[[ft_rmv]]$acc_time = calc_acc_time(match_cor_nd, "max")
  gps_ft_leison[[ft_rmv]]$acc_subj = calc_acc_subj(match_cor_nd, "max")
  print(mean(gps_ft_leison[[ft_rmv]]$acc_time))
}
  #permutation data
for (i in 1:15){
  ft_rmv = feat_names[i]
  print(ft_rmv)
  range = c(3:17)[-i]
  print(range)
  for (j in 1:perm_time_cut) {
    perm_part_times = 1
    print(paste("permuting ...", i,"..."))
    gps_df_perm$IID = sample(gps_df_perm$IID)
    perm_subj_seq = make_subj_seq(gps_df_perm, part_times = perm_part_times)
    perm_gps = make_feature_matrix(gps_df_perm,perm_subj_seq,range)
    perm_mat_1 = perm_gps$subj_mat_1
    perm_mat_2 = perm_gps$subj_mat_2
    perm_match_cor = calc_match_cor(perm_mat_1,perm_mat_2)
    perm_acc_time_cut[[ft_rmv]][[j]] = calc_acc_time(perm_match_cor)
    perm_acc_subj_cut[[ft_rmv]][[j]] = calc_acc_subj(perm_match_cor)
  }
  print(mean(unlist(perm_acc_time_cut[[ft_rmv]])))
}

acc_time_lesion_df$`None_Removed` = acc_time

Error in `$<-.data.frame`(`*tmp*`, None_Removed, value = c(0.48780487804878,  : 
  replacement has 1000 rows, data has 10

acc_time_lesion_perm_df <- as.data.frame(gather(acc_time_lesion_df, feature, acc_time, colnames(acc_time_lesion_df), factor_key=TRUE))

attributes are not identical across measure variables;
they will be dropped

subj_mins = gps_df_clean2 %>% group_by(IID) %>% dplyr::tally(name = "Days Collected")
conf_scatter_plot(subj_minmissing, gps_ft_leison$MinsMissing$acc_subj, "MinsMissing")

19. Misc Info

#sessionInfo()

 

---

A work by Cedric Huchuan Xia

hxia@upenn.edu